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INTRODUCTION 

This  report  contains  guidelines  for  using  the  circular  serial 
autocorrelation  filter  (CSAF)  program.  The  CSAF  program  is  a  later  version 
of  an  algorithm  called  SFILTR.  In  the  first  section  basic  things  the  user 
of  the  program  should  do  to  run  the  program  are  discussed.  The  requirements 
of  the  program  are  presented  in  this  first  section.  Interpretation  of 
the  program's  outputs  is  given  in  the  second  section.  In  the  next  section, 
things  which  the  user  should  not  do  are  indicated.  The  fourth  section  is 
a  three  part  discussion.  Mathematical  theory  related  to  the  filter  is  pre¬ 
sented  in  the  first  part;  then  the  program's  organization  in  terms  of  functions 
of  its  subroutines  and  main  program  is  discussed  in  the  second  part.  Follow¬ 
ing  this  part  are  listings  of  important  program  variable  definitions  and 
the  entire  program  along  with  a  flowchart  and  sample  run  streams  in  the  third 
part.  In  essence,  this  report  contains  the  most  detailed  description  of 
the  circular  serial  autocorrelation  filtering  tool  to  date. 
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1.  RUNNING  THE  CSAF  SOURCE  PROGRAM 


CSAF  is  a  numerical  algorithm  for  filtering  serial  data  utilizing  the 
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circularly  defined  autocorrelation  coefficient.  Real  data  to  be  filtered  must 
be  stored  in  an  existing  external  file  named  "INPUT"  to  be  list-directed  read 
by  the  source  program.  The  source  program  creates  an  estimate  of  the  input 
data  and  creates  a  corresponding  set  of  error  data.  In  total  before 
execution,  there  must  be  three  existing  files:  "INPUT",  "OUTPUT",  and 
"FHAT".  More  discussions  of  these  files  are  presented  in  Sections  2  and  3. 

CSAF  is  written  in  a  high  level  language,  FORTRAN.  A  FORTRAN  77 
compatible  compiler  is  needed  in  forming  an  executable  image.  The  main 
program  and  several  subroutines  reference  several  math  functions.  Thus,  the 
following  set  of  math  routines  also  should  be  available  for  double  precision 
applications:  sin,  cos,  arctan,  exp,  sqrt,  and  abs.  If  the  above 
requirements  are  met,  the  program  can  be  executed  on  good  personal  computers 
and  most  mainframes. 

It  is  assumed  that  a  set  of  parameters  related  to  the  input  signal's 
description,  the  desired  error  signal  statistical  test  constants,  and  the 
desired  estimated  signal  integration  method  are  known.  After  the  execution 
command  the  user  is  prompted  to  enter  six  parameters:  T,  NP,  GNT,  ZC,  NOT, 
and  NFLAG.  T  is  the  sampling  period  of  the  input  signal.  NP  is  the  number  of 
samples  of  the  input  signal.  GNT  is  the  Gaussian  noise  test  constant;  it  is 
the  critical  (alpha)  region  in  statistical  theory  used  in  testing  the  Gaussian 
nature  of  the  noise  signal.  ZC  is  the  critical  value  of  the  magnitude  of  the 
standard  normalized  Gaussian  random  variable;  it  is  used  in  testing  the 
randomness  of  the  error  signal.  NOT  is  a  logical  code  for  straight  line 


detrending  error  signals;  its  value  should  be  1  if  detrending  is  desired. 

NFLAG  is  a  logical  code  for  fast  Fourier  transforming  (FFT)  the  estimate 
signal  in  calculating  Fourier  coefficients:  its  value  should  be  1  if  an  FFT 
is  desired.  (Typical  values  for  GNT  and  ZC  are  0.05  and  1.96, 
respectively.)  After  entering  the  keyboard  input  parameters,  the  program 
should  produce  three  sets  of  outputs,  which  are  discussed  in  the  following 
section. 

2.  INTERPRETING  PROGRAM  OUTPUTS 

Outputs  from  the  source  program  are  written  to  the  terminal  screen,  into 
the  OUTPUT  data  file  and  into  the  FHAT  data  file.  The  terminal  displays  a  set 
of  pertinent  statements  which  provide  information  on  the  input,  estimate,  and 
error  signals.  These  statements  are  the  following: 

1.  A  statement  that  the  random  noise  is  non-Gaussian  at  the  alpha 
level,  if  the  test  for  normality  so  indicates; 

2.  An  estimate  of  the  signal-to-noi se  ratio  of  the  input  data; 

3.  The  program  determined  input  data's  average  value  (AO)  and 
fundamental  angular  frequency  (WO)  along  with  the  inputted  sampling 
period  (T),  number  of  data  points  (NP),  remaining  keyboard  input 
data;  and 

4.  The  even  and  odd  Fourier  coefficients  (As  and  Bs  respectively)  of 
the  estimated  signal . 


How  these  printed  outputs  are  calculated  is  discussed  in  Section  4.  Sample 
run  streams  are  given  on  pages  33  through  35. 

Real  output  data  are  stored  in  the  data  files.  The  OUTPUT  data  file 
contains  in  its  first  record  in  sequential  order  the  slope  of  the  straight 
line  detrend  (BETA),  curve  fit  error  (SE),  sampling  period  (T),  fundamental 
angular  frequency  (WO),  average  value  (AO)  and  optimum  cutoff  harmonic 
(L-1).  Remaining  records  of  OUTPUT  contain  the  even  and  odd  Fourier 
coefficients.  FUAT  contains  the  estimated  signal.  FHAT's  records  contain  a 
single  value  each,  and  are  sequentially  ordered  from  top  to  bottom,  the  top  of 
the  data  file  represent‘’‘ng  the  first  data  point  (initial  value  of  the  estimate 
signal).  In  some  applications,  either  mishandling  these  output  files  or  the 
input  file  presents  a  problem  in  running  the  program.  Reads  and  writes  from 
and  to  them  are  important  concerns  in  all  applications. 

3.  AVOIDING  MAJOR  PROBLEMS 

As  written,  all  external  files  are  direct  access  and  must  exist  before 
attempting  to  run  the  program.  Do  not  fail  to  create  the  INPUT,  OUTPUT  and 
FHAT  files.  The  program  assumes  real  values.  Use  of  complex  data  will  cause 
problems.  INPUT  must  contain  data  which  are  formatted  according  to  list- 
directed  requirements.  Do  not  place  more  than  one  entry  per  record  and 
neither  separate  entries  by  commas  nor  spaces.  If  a  subroutine  is  repeatedly 
called  in  a  redesigned  version  of  the  CSAF  program,  do  not  fail  to  reset  the 
associated  variables  (of  the  direct  access  files  affected)  to  the  first 
records  of  subroutine  files  called. 


In  addition  to  the  above  data  file  handling  pitfalls,  beware  of 
attempting  to  filter  long  data  records  greater  than  499  points,  especially  if 
it  is  known  that  the  data  have  been  poorly  sampled.  The  CPU  run  time 
increases  with  optimum  harmonic  content  of  the  data,  and  can  be  as  high  as 
2-1/2  hours  for  300  harmonics  of  2000  data  point  records  and  as  low  as  10 
seconds  for  20  harmonics  of  128  data  points.  Aliasing  cannot  be  completely 
removed  by  the  filter.  Aliasing  produces  cyclical  components  which  the 
statistical  method  assumes  are  not  present  and  the  program  cannot  remove 
through  the  subtraction  of  only  a  finite  number  of  harmonics. 

Other  concerns  are  the  following:  1)  the  number  of  input  data  points 
must  be  greater  than  75  for  asymptotic  characteristics  of  the  correlation 
statistics,  2)  disk  storage  for  CSAF  should  accommodate  24  blocks  (12288 
bits),  3)  the  Gaussian  noise  test  constant  should  not  be  zero  because  the  chi- 
square  tests  for  normality  will  not  converge,  and  4)  the  critical  z-value 
should  not  be  less  than  0.1  for  most  practical  problems.  More  explanations  of 
how  the  program  works  follow,  and  should  give  further  clarifications  of  why 
the  above  limits  are  necessary  for  good  program  execution. 

4.  UNDERSTANDING  THE  PROGRAM 

4.1  BASIC  THEORY 

Given  a  set  of  noise  contaminated  data,  the  objective  of  the  filtering 
routine  is  to  recover  the  continuous  desired  signal  without  thereby 
introducing  any  jump  discontinuities  in  the  recovered  signal  or  its 
derivatives.  Since  the  desired  signal  f(t)  often  contains  useful  or 


understandable  information,  it  is  often  referred  to  as  the  intelligence  signal 
in  the  sequel.  Because  of  the  possible  presence  of  random  noise  components  in 
discrete  recorded  data,  it  is  usually  not  possible  to  entirely  recover  the 
intelligence,  and  one  must  be  content  with  an  estimate  of  the  sampled  version 
of  this  signal,  fj^.  One  can  determine  an  intelligence  estimate  by  many 
methods.  In  this  work  the  estimate  signal  is  determined  indirectly  by 
construction  of  a  simulation  model  which  can  be  used  to  produce  a  sufficiently 
continuous  signal.  The  model  chosen  is  a  simple  truncated  Fourier  series 
whose  fundamental  frequency  and  coefficients  are  determined  from  the  input 
data.  Assuming  that  jump  discontinuities  of  the  noise  contaminated  continuous 
time  domain  signal  are  due  to  the  presence  of  random  noise  components,  the 
filtering  objective  can  be  met  by  utilizing  the  intelligence  estimate  to  get 
rid  of  most  random  noise  effects.  This  can  be  done  if  one  can  define  and 
isolate  the  noise  signal  or  its  estimate. 

Assuming  that  the  original  input  data  represents  an  evenly  sampled  signal 
whose  fundamental  period  is  the  data  length  (Np-1)  times  the  sampling  period 
(T).  i.e., 

Vq  -  (Np-1)  T  (1) 


0)^  -  2ii 
0  0 


(2) 


one  can  write  a  closed  expression  (model)  for  the  intelligence  estimate  for  L 
harmonics  at  every  data  point  k  as  the  following  summation: 
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(3) 


"  ^0  *  ^  [AjCOsCkiui^T)]  +  B^sin  (kitoj)], 

L  i”l 


1  i  L  ^  L. 


where  in  the  above  equation  k  e  {0,  N-1},  is  an  optimum  harmonic 

limit,  and  the  coefficients  are  given  for  i  t  {1,  L^}  through  the 

discrete  integrals 


1  f  ° 

r  !  I'(t)  dt 


0  0 


h  L 


e  ^ 

J  f(t)  COS  (iu^t)  dt,  and 


J  f(t)  sin  (luj^t)  dt. 


The  integrals  are  evaluated  using  a  fifth  degree  closed  Newton-codes  formula 
[1,  p.  142],  [2,  eq.  25.4.14],  Alternatively,  if  the  data  are  sampled  fast 
enough,  the  coefficients  can  be  obtained  through  fast  Fourier  transformation 
[3,  pp.  75-76].  Then,  one  can  define  the  error  signal,  having  produced  both 
intelligence  estimate  and  original  data  at  each  sample  point.  The  error 
signal  is  a  set  of  data  values  which  are  the  differences  between  the  input 
data  values  and  corresponding  intelligence  estimate  values,  i.e.. 


=kL„  •  , 

0  L, 


k  -  0,  ...,  N-1. 


Moreover,  one  can  define  a  family  of  error  signals  parametrized  by  the  maximum 
harmonic  content  (L)  of  the  intelligence  estimate.  Because  the  error  signal 
for  L  tends  to  become  more  random  as  L  increases,  the  estimate  of  the  random 
noise  component  is  defined  as 


k  -  0,  ...,  N-1. 


This  noise  estimate  can  be  removed,  and  the  resulting  signal  is  f^  —  as 
represented  by  equation  3  with  L  -  Lq-  In  addition  Xk  can  be  used  in 
calculation  of  a  signal-to-noise  ratio  estimate  for  the  contaminated  signal, 
This  figure  of  merit  is  given  by  the  following  expression; 


SNR  .  -20  log 


rms (x.  ) 

l\ 

rms  (f,^  ,  ) 
^0 


where  rms  means  “root  mean  square  value"  of  the  enclosed  quantity. 

Lq  in  the  above  paragraph  is  the  harmonic  at  which  the  curve-fit  error 
signal  ekL  (3S  a  function  of  L)  becomes  random*.  In  order  to  test  error 
signals  for  randomness,  their  pertinent  statistical  characteristics  are 
utilized.  First,  their  individual  means  (mL)  and  standard  deviations  (si_)  are 
determined;  then  each  rms  (e-^L)  is  found  by  the  formula 


rms  (e. ,  ) 


1  S  L  S  L. 


*  Lq  is  defined  to  be  less  than  N/2  to  prevent  violation  of  the  Nyquist 
sampling  rate. 
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And  from  these  rms  values,  curve-fit  error  norm  values  defined  by 


rms 


rms 


(11) 


are  determined  to  gauge  the  "closeness"  of  the  intelligence  estimates  to  the 
contaminated  signal  for  each  L.  Secondly,  the  circular  serial  autocorrelation 
coefficient  R  of  each  error  signal  and  the  Gaussian  nature  of  the  noise 
estimate  are  determined.  For  each  L-value  for  cyclically  defined  N**,  R  is 
given  by 


N-1 

«|l  -  ^N.L  V  ♦  -I, 


(12) 


where  tne  statistic  R  for  large  N  (N  >  75)  is  normally  distributed  having 
first  and  second  moments  respectively  given  by  the  following  [4]: 


(13) 


$2^  -  -  4S^^S2  +  43^83  +  $2^  -  2S4 

^  (N-1)  (N-2) 


(14) 


where  S^.(i-1,  ....  4)  is  the  ith  moment  about  the  origin  of  the  error 
signal.  Thus 


**  N  is  the  number  of  input  data  points  plus  an  extra  point  giving  f(NT)  - 


is  normal  with  zero  mean  and  unit  variance.  If  the  calculated  z-value 


magnitude  for  a  given  L  is  greater  than  1.96,  a  null  hypothesis  of  randomness 
can  be  confidently  rejected  at  the  5%  level.  Though  the  smaller  the  value  of 
jzj  the  more  confident  one  can  be  of  randomness,  in  most  practical  problems 
|Z|  is  not  zero;  but  may  be  much  less  than  1.96.  Various  critical  Izj-values 
can  be  used  for  various  levels  of  test.  1.96  is  simply  a  popular  one 
corresponding  to  a  5T  critical  region.  Having  |zl-value  for  each  L,  the  task 
is  simply  to  find  a  valley  turning  point  on  the  |z|  vs.  L  curve  below  the 
critical  z-value,  below  which  value  one  cannot  reject  a  hypothesis  of 
randomness . 

Now,  assuming  that  the  fist  random  error  has  been  determined  at  L  *  Lq, 
the  next  thing  to  do  is  to  determine  v/hether  the  hypothesis  that  the  resulting 
noise  is  Gaussian  is  true,  statistically.  This  can  be  done  by  performing  a 
chi-square  test  for  normality  (in  the  Gaussian  sense).  In  this  work  the 
dynamic  range  of  the  noise  estimate  (xj^)  is  divided  into  18  equally  spaced 
intervals  (groups);  then  the  theory  of  grouping  for  goodness  of  fit  is  applied 
[5,  pp.  85-90.]  The  cumulative  normal  distribution  for  standard  normalized 
noise  (zi<;)***  is  determined  using  a  fifth  degree  polynomial  —  no  relationship 


***  The  z.  variable  is  z. 


tc'" 


I 

II 


to  Newton-Cotes  —  approximation  [2,  eq.  26.2.17]  whicn  is  accurate  to  order 
10“^.  Theoretical  (t^)  and  observed  (0^)  frequencies  are  calculated.  The 
chi-square  statistic  is  determined  by 


^a,G-2  - 

X 


f-l  ^ 


where  a  is  the  level  of  significance,  G  is  the  number  of  groups  with  entries 
greater  than  four,  and  G-2  is  the  number  of  degrees  of  freedom.  One  can  now 
determine  the  area  under  a  chi-square  distribution  from  to  ®  using  the 
following  relation  for  v  degrees  of  freedom  [6,  p.  193]: 


♦  -  r  -1)  dX. 


2'"^r(\)/2)  2 

A 

However,  for  large  v  >  30  it  is  easier  to  use  the  inverse  formula 
approximation  for  Xa^  given  either  the  significance  level  or  (}>  [2, 
eq.  26.4.18];  and  to  use  Q-series  expansions  for  smaller  v  >  0  [2,  eqs.  26.4.4 
and  26.4.5],  i.e.. 


^  MWg  -  h  )J  ^  ] 


V  >  30 


_  2  j— — 2 

(2|v)  -  2  Q.(J7  )  +  2Z  (j7)  I 


V  odd  (19) 


Q(x  v) 


v-2 

2 

2^  fj?)  [1  +  I 

i-1 


-]  V  even 


Xa  should  be  greater  than  x  if  the  null  hypothesis  is  not  to  be  rejected.  The 
constant  h^  depends  on  the  value  of  Wq  [2,  eq.  26.4,15]  and  is  the  value  of 

the  standard  normal  random  variable  corresponding  to  the  a  critical  region. 

/ 

and  2|^J  have  usual  definitions  for  normal  random  variables: 


/I? 


Q,,  (u)  -  /  Z.(t)  dt. 


If  Q(x^|v)  is  greater  than  a,  the  hypothesis  of  normality  is  not  rejected. 
Otherwise,  it  is  rejected  and  the  non-Gaussian  message  is  printed  out  in  the 
computer  program. 

Finally,  trend  characteristics  of  the  error  signals  are  treated. 

Straight  line  trend  and  critical  movements  can  be  reduced.  Standard  equations 
for  the  abscissa  intercept  (alpha)  and  slope  (beta)  of  the  data  are  used  and 
are  not  repeated.  This  average  value  of  the  error  is  used  to  adjust  the  mean 
of  the  intelligence  estimate.  Cyclical  movement  is  minimized  by  subtracting 
from  the  error  signal  Fourier  harmonics  (by  design).  Thus,  the  error  should 
meet  the  hypotheses  of  Wald  and  Waltowitz  [4]  for  most  finite  data.  Having 
now  presented  the  essentials  of  the  theory  and  actual  use  of  the  filter 
program,  its  organization  and  examples  of  its  use  are  presented  next. 


4.2  PROGRAM  ORIENTATION 


The  circular  serial  autocorrelation  filtering  method  basically  involves 
determining  the  fundamental  frequency  of  noisy  data,  structuring  a  Fourier 
series  estimate  of  this  data  using  the  curve-fit  error  characteristics  in 
optimizing  the  Fourier  series  estimate  and  representing  the  "intelligence"  by 
a  set  of  Fourier  parameters.  These  basic  functions  are  performed  in  the  main 
program  of  CSAF.  In  the  first  part  of  the  program  the  number  of  data  points 
is  defined  for  circular  correlation,  and  the  fundamental  period  and 
fundamental  angular  frequency  are  computed.  Secondly,  the  Fourier  estimate  is 
computed  for  each  L  from  the  returned  values  of  Fourier  coefficients.  Next, 
the  curve-fit  error  is  calculated  for  each  intelligence  estimate.  Using  the 
correlation  statistic  the  turning  point  is  then  found  —  if  one  exists  below 
the  critical  z-value.  Finally,  the  set  of  Fourier  parameters  and  other 
information  describing  the  input  data  and  noise  estimate  are  printed  to  the 
screen.  In  performing  these  tasks  the  main  program  calls  several 
subroutines.  These  subroutines  are  discussed  in  the  following  paragraphs. 

To  begin  the  discussion,  there  are  two  options  as  indicated  by  the  values 
of  the  input  integers  NOT  and  NFLAG.  If  NOT  is  1,  the  DETRND  subroutine  is 
called.  This  routine  performs  least  squares,  straight  line  detrending  of  the 
error  signals  before  the  correlation  coefficient  of  the  noise  estimate  is 
calculated.  If  NFLAG  is  1,  the  FFT  subroutine  is  called  and  the  INTEGR 
subroutine  is  not  called.  The  FFT  subroutine  is  that  of  Bloomfield  [3]  with 
only  minor  adjustments;  its  function  is  the  calculation  of  the  Fourier 
coefficients.  If  NFLAG  is  not  1,  the  FFT  is  not  called  but  the  INTEGR  is 
called.  INTEGR  performs  the  ordinary  numerical  integration  of  the  waveform  to 


compute  the  Fourier  coefficients.  The  actual  integration,  however,  is 
performed  through  a  subsequent  call  to  the  QUAD  subroutine,  which  routine  is 
not  discussed.  With  NDT  and  NFLAG  inputted,  the  input  data  are  obtained 
through  a  call  to  the  DATAIN  subroutine.  This  routine's  function  is  simply  to 
read  INPUT,  set  the  last  data  point  equal  to  the  first  data  point  and  obtain 
input  data  statistics  from  a  call  to  SERCOR. 

The  most  frequently  called  subroutine  is  SERCOR.  This  subroutine 
performs  two  statistical  functions,  depending  on  the  value  of  its  fifth 
argument.  If  this  argument  is  2,  the  subroutine  simply  calculates  the  mean 
and  standard  deviation  of  the  signal  appearing  as  its  first  argument.  If  the 
fifth  argument  is  not  2,  the  subroutine  also  computes  the  circular  serial 
autocorrelation  of  the  signal  and  the  standard  normalization  of  this 
statistic.  There  are  three  other  subroutines  which  calculate  statistics  of 
signals.  The  GAUSS  subroutine  finds  the  complement  of  the  cumulative 
probability  [Q^]  of  a  normal  random  variable  using  a  polynomial 
approximation.  The  QSEXP  subroutine  finds  the  value  of  the  area  under  a  chi- 
square  probability  distribution  from  the  value  to  infinity  [Q(x^t'^)]  using 
the  proper  series  expansion  for  Q(x^l>>).  TNORM  performs  the  subgrouping  of 
the  noise  estimate  and  calculates  the  corresponding  chi-square  random 
variable.  Using  information  from  GAUSS  and  QSEXP,  TNORM  tests  the  hypothesis 
that  the  noise  estimate  is  random  at  the  GNT  level  of  significance.  Finally, 
the  TURN  subroutine  writes  data  to  OUTPUT  and  FHAT.  More  information  on  how 
the  program  functions  can  be  obtained  from  the  flowcharts  of  the  following 
section  (Figs.  1  through  3). 


4.3  PROGRAM  LISTINGS 


Following  is  a  list  of  important  variable  definitions  of  the  main 
program.  Also  included  in  the  alphabetical  list  are  important  variables  of 
the  TNORM  subroutine,  distinguished  from  the  main  program  variables  by 
asterisks.  This  list  is  not  complete  but  should  make  reading  the  source 
program  easier. 

A  Even  Fourier  coefficients,  an  array. 

ALP*  Critical  region  for  normality  test. 

ALPHA  Abscissa  intercept  of  straight  line  fit  of  error  signal. 

AVH  Average  value  of  input  data's  estimate. 

B  Odd  Fourier  coefficients,  an  array. 

BETA  Slope  of  straight  line  fit  of  error  signal. 

CHIS*  Chi-square  statistics  of  noise  estimate. 

CNP  Cumulative  normal  probability  of  standard  normalized  noise 
estimate,  an  array. 

CX*  Intermediate  value  of  chi-square  statistics  of  noise  estimate,  an 
array. 

DOF*  Degrees  of  freedom  for  noise  estimate  chi-square  statistic. 

DX*  Width  of  groups  of  making  up  noise  estimate  signal. 

DDX*  Discrete  boundary  of  DX  groups,  an  array. 

E  Error  norm,  an  array. 

EAVG  Average  value  of  error  signal. 

ERR  Error  signal,  an  array. 

ESIG  Standard  deviation  of  error  signal. 
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FN 


Input  data,  an  array. 

FH  Intelligence  (Fourier)  estimate,  an  array. 

FQ*  Cumulative  chi-square  probability  of  noise  estimate  chi-square 
statistic. 

GMT  Gaussian  noise  test  constant. 

LMAX  Maximum  number  of  harmonics  calculated. 

N  Number  of  points  in  circular  data. 

ND*  Number  of  sections  into  which  the  range  of  noise  signal  is 
divided. 

NOT  Logical  code  for  straight  line  detrending  of  error  signal. 

NFLAG  Logical  code  for  performing  FFT. 

NP  Number  of  input  data  points. 

NR  Logical  code  for  calculation  of  correlation  coefficient  of  data. 

OF*  Observed  frequency  of  values  of  noise  estimate  within  OX  limits, 

an  array. 

0*  Cumulative  normal  probability  of  standard  normalized  noise. 

SAH  Average  value  of  input  data's  estimate. 

T  Sampling  period  of  input  data. 

TF  Theoretical  frequency  of  values  of  noise  estimate  within  DX 
limits,  an  array. 

TRF  Theoretical  relative  frequency  of  noise  estimate  values,  an  array. 
WO  Fundamental  angular  frequency  of  input  data. 

X  Detrended  error  signal,  an  array. 

X*  Noise  estimate  signal,  an  array. 
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XI 


Imaginary  part  of  input  data  (zero)  and  odd  Fourier  coefficient  in 
FFT  return,  an  array. 

XR  Real  part  of  input  data  and  even  Fourier  coefficient  in  FFT 
return,  an  array. 

Z  Standard  normalized  correlation  coefficient. 

Z*  Standard  normalized  noise  estimate  values,  an  array. 

ZC  Critical  magnitude  of  Z-value. 

ZX  Z-values  of  family  of  error  signals,  an  array. 

A  copy  of  the  entire  source  program  is  given  on  pages  19  through  26. 
Flowcharts  for  the  main  program  and  the  SERCOR  and  TNORM  subroutines  are  given 
on  pages  27  through  32.  Sample  run  streams  are  given  on  pages  33  through  35. 

5.  CLOSING  REMARKS 

We  have  attempted  to  answer  some  fundamental  questions  concerning  the 
computer  program  used  to  low-pass  filter  data  based  on  the  serial  correlation 
coefficient  of  curve-fit  error.  Though  we  have  addressed  the  program's 
inputs,  outputs  and  principles  of  operation,  there  can  never  be  a  complete 
description  of  the  intricacies  of  efficient  utilization  of  the  program.  This 

happens  after  using  the  program  a  number  of  times  and  applying  variations  of 

the  basic  theory  on  which  the  algorithm  is  based  to  practical  problems.  We 
have  done  this  using  various  experimental  and  simulated  data  and  we  have  not 
been  able  to  find  a  case  where  the  program  has  not  done  what  it  is  designed  to 
do.  We  welcome  all  comments  on  the  performance  of  the  program  from  all  users. 
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CSAF  SOURCE  PROGRAM 


C  THIS  PROGRAM  USES  SERIAL  CORRELATOR  METHOD  TO  FILTER  INPUT 
C  DATA.  FOURIER  COEFFICIENTS  OF  FILTERED  DATA  ARE  FOUND  UP 
.C  .  TO  AN  OPTIMUM  NUMBER  OF  HARMONICS.  THE  SIGNAL  TO  NOISE  RATIO 

C  *  OF  THE  ORIGINAL  INPUT  DATA  IS  APPROXIMATED  AND  AN  ERROR 

C  ESTIMATION  OF  THIS  METHOD  IS  COMPUTED. 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(250),B(250),X(500),FN(500),FH(500) 

DIMENSION  XR(500),XI(500) 

DIMENSION  ERR(500),ZX(250),E(500) 

C  THIS  SECTION  READS  THE  SET  OF  INPUT  PARAMETERS 
C  THE  INPUT  PARAMETERS  REQUIRED  ARE: 

C  T  :  SAMPLING  PERIOD. 

C  NP  :  NUMBER  OF  SAMPLED  POINTS  TO  BE  USED. 

C  GNT  ;  GAUSSIAN  NOISE  TEST  CONSTANT. 

C  ZC  :  CRITICAL  VALUE  OF  Z. 

C  NDT  :  DATA  DETRENDING  CODE.  NDT=1  IF  DETREND  IS  TO 

C  BE  USED. 

C  NFLAG  :  FFT  CODE.  IF  FFT  IS  TO  BE  USED,  NFLAG=1. 

PI=4.0*DATAN(1.D0) 

ZX(l)=O.DO 

WRITE(6,245) 

10  READ(5,*)T,NP,GNT,ZC,NDT,NFLAG 
VRITE ( 6 , 2  60 ) GNT , ZC , NFLAG , NDT 
N=NP+1 

IF(N.LE.75)WRITE(6,200) 

IF(N.LE.75)G0  TO  300 
CALL  DATAIN(N,FN,AV,SA) 

IF(NFLAG.EQ.l)  GO  TO  50 
17  DO  20  1=1, N 
20  ERR(I)=FN(I) 

50  W0-2.D0*PI/((N-1)*T) 

LMAX=(N-l)/2 

LTN=0 

IF(NFLAG.EQ.l)  GO  TO  30 
C  NP=.5DO+2.DO*PI/(T*WO) 

C  INT=N/NP 

GO  TO  40 

C  THIS  SECTION  DETERMINES  THE  FOURIER  COEFFICIENT  OF  THE 
C  INPUT  DATA  USING  FAST  FOURIER  TRANSFORM  (FFT). 

30  DO  170  I-1,N 
XI(I)-O.DO 
170  XR(I)-FN(I) 

CALL  FFT(XR,XI,N-1,0) 

AO-XR(l) 

DO  180  I-1,LMAX 
A(I)-2.D0*XR(I+1) 

B(I)— 2.D0*XI(I+1) 

180  CONTINUE 

C  THIS  SECTION  RECONSTRUCTS  THE  ORIGINAL  INPUT  DATA  USING 
C  FOURIER  COEFFICIENT  APPROXIMATION  AND  ESTIMATES  THE 
C  ERROR  VALUES.  IF  FFT  IS  NOT  USED,  A  NUMERICAL  INTEGRATION 
C  METHOD  IS  USED  TO  DETERMINE  THE  FOURIER  COEFFICIENTS. 

40  DO  190  L-1,LMAX 
NR-0 

IF(NFLAG.EQ.l)  GO  TO  60 

CALL  INTEGR(T,N,UO,AO,A,B,L,ERR,NP) 

60  DO  100  K-1,N 


non 


•  AUG=I*W0*(K-1)*T 

70  FE=FS+A(I)*DCOS(AUG)-rB(I)*DSIN(AUG) 

FH(K)=.FE+A0 
ERR(K)=FN(K)-FH(K) 

100  X(K)=ERK(K) 

THIS  SECTION  USES  THE  SERIAL  CORRELATOR  METHOD  TO  FILTER 
THE  INPUT  DATA  AND  ESTIMATES  THE  S/N  RATIO.  FOURIER 
COEFFICIENTS  OF  THE  FILTERED  DATA  ARE  OUTPUTTED. 

CALL  SERC0R(FH,N,SAH,AVH,2,AZ) 

CALL  SERCORC ERR , N , ESIG , EAVG , 2 ,  EAZ) 

I F ( NDT . EQ . 1 ) CALL  DETRND ( X , N- 1 , ALPHA , BETA ) 

120  CALL  SERC0R(X,N,SIG,AVG,NR,2) 

E(L)=DSQRT((ESIG**2+EAVG**2)/(SAH**2+AVH**2)) 

ZX(L-t-l)=Z 

IF(L.LT.3)G0  TO  125 

IF((ZX(L).GE.ZX(L-1)).0R.(ZX(L).GE.ZX(L+1)))G0  TO  125 
IF(ZX(L).GE.ZC)GO  TO  125 
CALL  TNORM(GNT,SIG,AVG,N,X,IFUG) 

VRITE(6,*)'  SNR  <=>  ' , -2 .D1*DL0G10(E(L) ) , '  DB' 
IF(IFLAG.EQ.1)URITE(6,290)GNT*100. 

CALL  TURN(N,L,FH,SEE,T,VO,AO,A,B,ZX,BETA,ALPHA,LTN) 

GO  TO  195 
125  CONTINUE 
SEE=E(L) 

190  CONTINUE 
195  CONTINUE 

VRITE(6,270)SEE 
IF(LTN.NE.1)VRITE(6,275) 

WRITE (6, 207) 

URITE(6,210)A0,W0,T,N-1 
VRITE(6,215) 

DO  160  1=1, L-1 
160  VRITE(6, 220)1, A(I),B(I) 

VRITE(6,250) 

300  CONTINUE 
310  STOP 

200  FORMATC'O' ,'N  IS  LESS  THAN  75.  PROGRAM  ABORTS.') 

207  FORMATC'O' ,/,6X,' AO' ,11X, 'WO' ,12X,'T',11X,'N' ,/) 

210  FORMAT('  ' ,3D13.6,I6,/) 

215  FORMATC'O' , 'HARMONIC' ,10X, 'A' ,15X, 'B' ,/) 

220  FORMATC'  ' ,3X,I3,5X,D13.6,5X,D13.6) 

230  FORMATC'  ' ,D13.6,12X,D13.6) 

245  FORMATC IX, 'ENTER  T,NP,GNT,ZC,D-CODE,FFT-CODE' ,/) 

250  FORMATC//) 

260  FORMATC'O' , 'GAUSSIAN  NOISE  TEST  CONSTANT  -',F6.4,/, 

&1X, 'CRITICAL  Z-VALUE  =' ,F7 .4,/ ,1X, ' FFT-CODE  -' , 12 , / , IX, 'D-CODE 
&I2,/) 

270  FORMATC 'O', 'RMS  CURVE-FIT  ERROR  =',E12.5,/) 

275  FORMATC'O' , 'TURNING  POINT  BELOW  CRITICAL  Z-VALUE  WAS  NOT  FOUND. 
&//) 

290  FORMATC'O' ,' NOISE  IS  NON-GAUSSIAN  AT  THE',F4.1,'  Z  LEVEL.') 

END 

C . 

c . 

SUBROUTINE  DETRND C Z, L, ALPHA, BETA) 

IMPLICIT  REAL*8CA-H,0-Z) 

DIMENSION  ZC0:499) 

REAL*8  A,B,C1,C2,DELTA,T1,T2,S1,S2,TS 
REAL*8  XC0:499),YC0:499) 

10  DO  90  K-1,2 
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S2=0.D0 

S1=0.D0 

T1=0.D0 

T2=0.D0 

TS=O.DO 

DO  20  1=0, L 

IF(K.EQ.l)GO  TO  15 

Y(I)=-C2*X(I)+C1*(Y(I)-A) 

X(I)=C1*X(I)+C2*(Y(I)-A) 

CONTINUE 

IF(K.EQ.1)X(I)=DFLOAT(I+0) 

IF(K.EQ.1)Y(I)=Z(I) 

T1=T1+X(I) 

T2=T2+X(I)**2 

TS=TS+Y(I)*X(I) 

S1=S1+Y(I) 

S2=S2+Y(I)**2 

DELTA=(L+1)*T2-T1**2 

A=(T2*S1-T1*TS)/DELTA 

B=((L+1)*TS-T1*S1)/DELTA 

VM=S1/(L+1) 

SIG=DSQRT((S2-S1**2/(L+1))/L) 
XTY=(L+1)*A-S1+B*T1 
IF(K.NE.1)G0  TO  80 
C1=DC0S(DATAN(B)) 
C2=DSIN(DATAN(B)) 

ALPHA=A 
BETA=B 
GO  TO  90 
DO  50  1=0, L 
2(I)=Y(I) 

CONTINUE 

CONTINUE 

RETURN 

END 


C. . . . 
C.  •  •  • 


SUBROUTINE  TNORM ( ALP , XSIG , XM , N , X , IFLAG ) 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  DDX( 25 ) , Z ( 25 ) , X(N) , CNP( 25 ) , CX( 25 ) 
DIMENSION  TF(25),TRF(25) 

INTEGER  0F(25) 

XMAX=X(1) 

XMIN=X(1) 

ND=20 

IFLAG=0 

CHIS=0.D0 

CNP(1)=0.D0 

CNP(ND+1)=1.D0 

DDX(l)— 10.D36 

DDX(ND+1)=10.D36 

DO  30  1=1, N 

IF(X(I).GE.XMAX)XMAX=X(I) 

IF(X(I).LE.XMIN)XMIN-X(I) 

CONTINUE 

DX=(XMAX-XMIN)/(ND-2) 

DO  100  I=1,ND-1 
DDX(I+l)=DX*(I-ND/2)+XM 
Z(I)-(DDX(I+1)-XM)/XSIG 
A-Z(I) 

CALL  GAUSS(A,Q,AZ,AS2PI) 

CNP(I+1)-1.D0-Q 

CONTINUE 

DO  120  J=1,ND 


rW.w^ 


>_V,- 


no 


TRF(J)=CNP(J+1)-CNP(J) 

TF(J)=TRF(J)*N 
•  0F(J)=0 

DO  113  1=1, N 

IF((X(I).GT.DDX(J)).AND.(X(I).LT.DDX(J+1)))0F(J)= 

&1+0F(J) 

118  CONTINUE 
120  CONTINUE 
IC=0 
IB=0 

130  DO  200  J=1,ND 
132  CONTINUE 

IF((OF(J).GE.5).OR.(J.EQ.ND))GO  TO  195 
0F(J+1)=0F(J)+0F(J+1) 

TF(J+1)=TF(J)+TF(J+1) 

IC=IC+1 

140  DO  180  K=1,ND+1-IC-J 
0F(J+K-1)=0F(K+J) 

TF(J+K-1)=TF(K+J) 

180  CONTINUE 

IF((J.GT.ND-IC).AND. ((IB+IC).LE.ND))IB=IB+1 
IF(J.GT.ND-IC)GO  TO  210 
GO  TO  132 
195  1B=IB+1 
200  CONTINUE 
210  D0F=IB-3 
220  DO  240  J=1,IB 

CX(J)=((0F(J)-TF(J))**2)/TF(J) 

CHIS-CHIS+CX(J) 

240  CONTINUE 

CALL  QSEXP(DOF,CHIS,FQ) 

AT=1. DO-ALP 
IF(FQ.GE.AT)IFLAG=1 
300  RETURN 
END 


SUBROUTINE  QSEXP(XNU,XO,FQ) 
IMPLICIT  REAL*8(A-e,0-Z) 

10  X=DSQRT(XO) 

CALL  GAUSS(X,Q,Z,S2PI) 

NU=XNUfl.D-9 

JEFF.INT(XNU/2.D0) 

MUTT=INT( XNU-2 . D0*JEFF+ . 5D0 ) 

QSU=0.D0 

ID-1 

IF(MUTT.EQ.0)G0  TO  140 
100  DO  120  I-l,(NU-l)/2 
ID-<2*I-1)*ID 

120  QSU-QSU+(X**(2*I-1))/ID 
QXN-2 . D0*Q+2 . D0*Z*QSU 
GO  TO  160 

140  DO  150  I=2,(NU-2)/2 
ID-2*I*ID 

150  QSU-QSU-^(X**(2*I))/ID 
QXN-S2PI*Z*(1.D0+QSU) 

160  FQ-l.DO-QXN 
200  RETURN 
END 


SUBROUTINE  GAUSS(X,Q,Z,S2PI) 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  B(5) 
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DATA  (3(1), 1=1, 5)  /.31938153CDO, -.35656378200,1.78147793700 
£.-1.82125597800, 1.33027442900/ 

PI  =  4.D0*DAT.AN(1.D0) 

S2PI=DSQRT(2.0Q*?I) 

X0=X-**2 

Z=(DEXP(-XO/2.O0))/S2PI 
XPP=. 231641900 
T=1.D0/(1.D0+X?P*X) 

SUM=0.D0 
50  DO  75  1=1,5 
75  SUM=SUM-3(I)’^(T^*(I)) 

80  Q=Z*SUM 
150  RETURN 
END 


SUBROUTINE  DATAIN(N,F,AVG,SIG) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  F(500) 

OPEN  ( 30, FILE=' INPUT. DAT' , STATUS=' OLD' ) 
DO  10  1=1,126 
10  READ(30,*)  F(I) 

CALL  SERCOR(F,N,SIG,AVG,2,O.DO) 
F(N)=F(1) 

80  RETURN 
END 


SUBROUTINE  INTEGR(T , N , VO , AO , A , B , NFQ , F ,  INT ) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  F(N) ,A(NFQ) ,B(NFQ),D(500) ,E(500) ,U(500) ,V(500) 

K=NFQ 

U(1)=0.D0 

V(1)=0.D0 

D(1)=F(1) 

E(1)=0.D0 
50  WK=(K)*W0*T 
DO  120  L=2,INT+1 
AUG=WK*(L-1) 

D(L)=F(L)*DCOS(AUG) 

90  E(L)=F(L)'^DSIN(AUG) 

120  CONTINUE 

CALL  QUAD(IC,INT,T,D,SIA) 

U(INT+1)=SIA 

CALL  QUAD(IC,INT,T,E,SIA) 

V(INT+1)=SIA 

A(K)=2.D0*U(INT-hl)/(T*INT) 

B  (  K ) = 2 .  DO*  V  ( INT-e  1 )  /  (  T*INT ) 

IF(K.GT.1)G0  TO  160 
A0=.5D0*(F(1)+F(INT+1))/INT 
CALL  QUAD(IC,INT,T,F,AO) 

A0=A0/(INT*T) 

160  CONTINUE 
210  RETURN 
END 


SUBROUTINE  QUAD(IC,INT,T,Q,SIA) 
IMPLICIT  REAL*8(A-H,0-Z) 
DIMENSION  Q(500) 

SIA=0.D0 

IC=0 

DO  130  I=l,INT-4,4 


SIA=SIA+2.0D0*T*(7.0D0*Q(I)+32.D0*Q(I+l)+12.D0*Q(I+2) 
&+32.D0*Q(I+3)+7.D0*Q(I+4))/45.D0 
•  .130  CONTINUE 

IH=INT-4*IC-1 

IG=1-4*IC 

IF(IH.EQ.3)SIA=SIA+3.D0*T*(Q(IG)-i-3.D0*Q(IG+l)+3.D0*Q(IG+2) 

5,+Q(IG+3))/8.DO 

IF(IH,EQ.2)SIA=SIA+T*(Q(IG)+4.D0*Q(IG+l)+Q(IG+2))/3.D0 
IF(IH.EQ.1)SIA=SIA+T*(Q(IG)+Q(IG+1))/2.D0 
160  RETURN 
END 

C . 

C . 

SUBROUTINE  SERCOR ( X , N , SIG , AVG , NR , Z) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  X(N),S(4) 

R=O.DO 
DO  10  1=1,4 
10  S(I)=O.DO 
DO  40  1=1,4 
DO  40  J=1,N 
40  S(I)=S(I)+X(J)**I 
AVG=S(1)/N 

SIG=DSQRT ( ( S ( 2 ) -N*AVG**2 ) / (N-1 ) ) 

IF(NR.EQ.2)G0  TO  100 
EXR=(S(1)**2-S(2))/(N-1) 

VAR=(S(2)**2-S(4))/(N-1)-(S(1)**2-S(2))**2/(N-1)**2 

VAR=VAR+(S(l)**4-4*S(2)*S(l)**2+4*S{l)*S(3))/((N-l)*(N-2)) 

VAR=VAR+(S(2)**2-2*S(4))/((N-l)*(N-2)) 

DO  70  1=1, N-1 
70  R=R+X(I)*X(I+1) 

R=R+X(1)*X(N) 

Z=DABS(R-EXR)/DSQRT(VAR) 

IF(Z.GE.25.D-4)NR=1 
100  RETURN 
END 


C . 

C . 

SUBROUTINE  FFT(XR,XI,N,INV) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  XR(N) , XI (N) , UR ( 15 ) , UI ( 15) 

LOGICAL  FIRST 

DATA  FIRST  /.TRUE./ 

IF(.NOT.FIRST)GO  TO  120 
UR(1)=0.D0 
UI(1)=1.D0 
DO  110  1=2,15 

UR(I)=DSQRT((1.D0+UR(I-1))/2.D0) 

110  UI(I)=UI(I-1)/(2.D0*UR(I)) 

FIRST=. FALSE. 

120  IF(N.GT.O  .AND.  N. LE. 2**16)GO  TO  130 
INV=-1 
RETURN 
130  N0=1 
II-O 

140  NO=NO+NO 
II  =  Il4-l 

IF(NO.LT.N)GO  TO  140 

Il-NO/2 

13-1 

lO-II 

DO  260  14-1,11 
DO  250  K=1,I1 


n  o 


UI=O.DO 

KK=K-1 

DO  230  1=1,10 
IF(KK.EQ.O)GO  TO  240 
IF(M0D(KK,2).EQ.0)G0  TO  230 
JO=IO-I 

WS=UR*UR(JO)-WI*UI(JO) 

WI=UR*UI(JO)+UI*UR(JO) 

VR=WS 

230  KK=KK/2 

240  IF(INV.EQ.O)VI=-yi 
L=K 

DO  250  J=1,I3 
L1=L+I1 

ZR=XR(L)+XR(L1) 

ZI-XI(L)+XI(L1) 

Z=VR*(XR(L)-XR(L1))-UI*(XI(L)-XI(L1)) 
XI(L1)=WR*(XI(L)-XI(L1))+VI*(XR(L)-XR(L1)) 
XR(L1)=2 
XR(L)=ZR 
XI(L)=ZI 
250  L=L1+I1 
10=10-1 
13=13+13 
260  11=11/2 
UM=1.D0 

IF( INV . EQ . 0) UM=1 . DO/DFLOAT(NO) 

DO  310  J=1,N0 

K=0 

J1=J-1 

DO  320  1=1,11 
K=2*K+M0D(J1,2) 

320  Jl=Jl/2 
K-K+1 

IF(K.LT.J)GO  TO  310 
ZR=XR(J) 

ZI=XI(J) 

XR(J)=XR(K)*UM 
XI(J)=XI(K)*(JM 
XR(K)=ZR*UM 
XI(K)=ZI*UM 
310  CONTINUE 
400  RETURN 
END 


•  •  •  •  « 

SUBROUTINE  TURN(N,L,FH,SE,T,WO,AO,A,B,ZX,BETA,ALPHA,LTN) 
IMPLICIT  REAL*8(A-a,0-Z) 

DIMENSION  FH<N) ,A(L) ,B(L) ,ZX(L+1) 

OPEN  (20,FILE='0UTPUT.DAT' ,STATUS='OLD' ) 

OPEN  (21, FILE-' FHAT.DAT' , STATUS-'OLD' ) 

LTN-1 

AO=AO+ALPHA 

WRITE(20,*)  BETA,SE,T,W0,A0,L-1 
DO  80  1=1, L-1 
80  VRITE(20,*)A(I),B(I) 

VRITE(21,*)(FH(I),I=1,N) 

200  RETURN 
END 
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V-r.V^ 


Figure  1 
CSAF 


IF  " 
NFLAG=1 


/  DO  (— 

L=1 ,LMAX  ^  C 


(A,B,L,X,NP) 


CALL  SERCOR 
(FH,N,SA,AV, 
2,AZj _ 


GALL  SERCOR 

(err;n,esig, 

EAUG,2~EAZ) 


Figure  1  (Continued) 


I-l  ,N 


XI(I)=0 

XR(I)=FN(I) 


Call  FFT(XR, 
XI,N-1 .0) 


A(I)=2*XR(I) 

B(I)=2*XI(I) 


'DO 

1=1. L 


FH(K)=FE+A0 
ERR(K)=FN(k) 
-FH  k) 
X(k)=  ERR(K) 


AUG=I*W0*T 

*(k-l) 


■FE=FE+A(I)* 
COS(AUG)  ! 
+B(I)*SIN(AUG] 


FIGURE  2 

Subroutine  SERCOR 


S(2)**2-S*S(A) 

(N-i)=HN-:) 


FIGURE  3 


SAMPLE  RUN  STREAMS 


ERD131S  RUN  CSAF 

ENTER  T,NP,GNT,2C,D-C0DE,FFT-C0DE 


1.949147E-10,128, .05,1.96,1,0 
GAUSSIAN  NOISE  TEST  CONSTANT  -0.0500 
CRITICAL  Z-VALUE  -  1.9600 
FFT-CODE  -  0 
D-CODE  -  1 


SNR  <=>  7.988648453179364  D3 
NOISE  IS  NON-GAUSSIAN  AT  THE  5.0  Z  LEVEL. 
RMS  CURVE-FIT  ERROR  =  0.40421E+00 


AO  VO  T 

0.135005D+C6  0. 251840D-f.09  0.194915D-09 


HARMONIC 


FORTRAN  STOP 


B 


1 

-0.196020D+07 

-0.4374640+06 

2 

-0.167392D+07 

-0.1309780+07 

3 

-0.144792D+07 

-0.1148140+07 

4 

-0.217603D+07 

-0.2382420+07 

5 

-0. 1739920+07 

-0.4022660+07 

6 

-0. 4004720+06 

-0.6816200+07 

7 

0.3643140+07 

-0.9515000+07 

8 

0.1127350+08 

-0.2988690+07 

9 

0.8007900+07 

0.1900040+07 

10 

0.5269640+07 

0. 3159660+07 

11 

0.3967800+07 

0.2609170+07 

12 

0.3098130+07 

0.2530830+07 

13 

0.3001510+07 

0.2352670+07 

14 

0.2375610*07 

0. 2068480+07 

15 

0.2134470+07 

0. 1928780+07 

16 

0.2127990+07 

0.1861680+07 

17 

0.1890280+07 

0.1669670+07 

13 

0.1991060+07 

0.1673980+07 

19 

0.2014040+07 

0.1733480+07 

20 

0.1895100*07 

0.1783780+07 

21 

0.2065990+07 

0.1990340+07 

22 

0.1967440+07 

0.2638100+07 

23 

0.1349000+07 

0.3462950+07 

24 

-0.6054330*06 

0.3585830+07 

25 

-0.1470680+07 

0.2161880+07 
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ERD131$  RUN  CSAF 

ENTER  T , NP , GNT , ZC , D-CODE , FFT-CODE 

1. ,128, .05,1.96,1,0 

GAUSSIAN  NOISE  TEST  CONSTANT  =0.0500 

CRITICAL  Z-VALUE  -  1.9600 

FFT-CODE  =  0 

D-CODE  =  1 

R.MS  CURVE-FIT  ERROR  =  0.86745E+00 

TURNING  POINT  BELOV  CRITICAL  Z-VALUE  WAS  NOT  FOUND. 


AO 

VO 

T 

-0.156881D- 

-01  0.490874D-01 

O.lOOOOOD+01 

HARMONIC 

A 

B 

1 

-0.2244000-01 

-0.1534700-01 

2 

-0.195964D-01 

-0.1771220-01 

3 

-0.108471D-01 

-0.2673600-01 

4 

-0.316106D-02 

-0.1929210-01 

5 

-0. 5324860-03 

-0.1711390-01 

6 

0.1739190-02 

-0.1405270-01 

7 

0.4905720-02 

-0.1244230-01 

8 

0.633246D-02 

-0.841825D-02 

9 

0.6395620-02 

-0.6618880-02 

10 

0.7450680-02 

-0.5453580-02 

11 

0.6673890-02 

-0. 3086890-02 

12 

0.6494230-02 

-0.1977640-02 

13 

0.4393870-02 

-0. 1284270-02 

14 

0.5062330-02 

0.4765190-04 

15 

0.4379110-02 

0.3408980-03 

16 

0.2663660-02 

0.4421730-03 

17 

0.2544130-02 

0.1160150-02 

18 

0.1032610-02 

0.1069050-02 

19 

0.2555650-03 

0.1216040-02 

20 

-0.6687290-03 

0.8520490-03 

21 

-0.1875930-02 

0.8301230-03 

22 

-0.3375330-02 

-0.1390400-02 

23 

-0.4243100-02 

-0.2272550-02 

24 

-0.5038190-02 

-0.3928230-02 

25 

-0.3495500-02 

-0.5707550-02 

26 

-0.3148270-02 

-0.4675190-02 

27 

-0.5807680-02 

-0.5027370-02 
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28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 


-0.100394D-01 
-0. 4949630-02 
0.319914D-02 
0. 1009480-01 
0.1534150-01 
0.2619580-01 
0.2517880-01 
0.2466310-01 
0.1662330-01 
0.6991550-02 
0.6746090-03 
0.1007680-02 
0.3734250-02 
0.5318470-02 
■  0.6015530-02 
0.6235310-02 
0.6256360-02 
0.6033890-02 
0.6495170-02 
0.6363160-02 
0.6838050-02 
0.5998090-02 
0.5669250-02 
0.4980970-02 
0.3556200-02 
0.2828230-02 
0.2044910-02 
0. 8113760-03 
0. 9827430-04 
-0.2970290-03 
-0.5801230-04 
0.1964780-02 
0.5034450-02 
0.9953900-02 
0.1472010-01 
0.1211930-01 
0.1627750-01 


-0. 1214690-01 
-0. 1995980-01 
-0. 2362300-01 
-0.2115100-01 
-0.2207440-01 
-0. 6700050-02 
-0.2384650-02 
0.3064430-02 
0. 1177850-01 
0. 1427410-01 
0.7573960-02 
0.8485100-03 
-0.3431990-03 
-0.5334710-03 
-0.6694130-03 
0.8843200-03 
0.6282630-03 
0.1162710-02 
0.1030420-02 
0.1295420-02 
0.203499D-02 
0.2761880-02 
0.3563140-02 
0.3351970-02 
0.3278400-02 
0.2887590-02 
0.2434240-02 
0.1234260-02 
-0.1860410-03 
-0.1900780-02 
-0.4721910-02 
-0. 7154460-02 
-0.9369450-02 
-0. 8784760-02 
-0.4411180-02 
-0.1174080-02 
-0.1424280-16 


FORTRAN  STOP 
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